function initialize ( datafile )
%% data handling for the optimization
% Read data file 'case$x.m'
name = mfilename('fullpath'); 
procpath = name(1:max([strfind(name,'/') strfind(name,'\')]));
global dataDir dataID
dataDir = strcat(procpath, 'case\'); dataID = strtok(datafile, '.');
addpath(dataDir);

global bus gen branch genpara sbspara
% Execute '.m' files
eval('global NTswitch');            % To remove warnings
run(strcat(dataDir, datafile));
% Read option file for cases
caseOption(datafile);
rmpath(dataDir);

global nB nL nG   iLD iG iG_cut iGB_cut   signG uGez
global TJ TJsum Xd1   idx_gx idx_gy
global Gld_cut Bgen_cut Bld_cut   Pgen_cut Pld_cut
global omgs omega delta Pm EXdTJ EXd
% Elimilate --- the off-line bus, generator and branch
bus(bus(:,2)>=4, :) = [];                  % find off-line buses[logical] and clear. 
ig = gen(:,6)==0;                          % find off-line generators
gen(ig, :) = []; genpara(ig, :) = [];      % clear off-line generators
branch(branch(:,11)==0, :)=[];             % clear unconnected branches

% Standardization --- power per unit value
bus(:, 4)=bus(:,4)/baseMVA;
bus(:, 5)=bus(:,5)/baseMVA;
bus(:, 6)=bus(:,6)/baseMVA; 
bus(:, 7)=bus(:,7)/baseMVA; 
bus(:, 9)=bus(:,9)*pi/180;         % Vangle of bus: degree to radius
gen(:, 2)=gen(:,2)/baseMVA;
gen(:, 3)=gen(:,3)/baseMVA;
gen(:, 8:11)=gen(:,8:11)/baseMVA;

branch(:,6)=branch(:,6)/baseMVA;
branch(:,7)=branch(:,7)/baseMVA;
branch(:,8)=branch(:,8)/baseMVA;
branch(:,10)=branch(:,10)*pi/180;  % ?

%% some parameters of the given system
nB = size(bus,1);                  % buses(system nodes)
nL = size(branch,1);               % branches
nG = size(gen,1);                  % generators

iLD = find(bus(:,3)==1);           % cuttable load bus
iG = gen(:,1);                     % node number of generators
signG = gen(:,6);%gen(:, 7);
uGez = find(gen(:,7)==0);
iG_cut = find(gen(:,6)~=0);%find(gen(:,7)~=0);        % gen number of cuttable generator
iGB_cut = gen(iG_cut, 1);          % node number of cuttable generator

idx_gx = iG; idx_gy = iG + nB;
TJ = genpara(:,7) * 2; TJsum = sum(TJ);
Xd1 = genpara(:,6);                % Xd'
omgs = 2*pi*sbspara(3);            % omega_s
%% Shunt admittance of cuttable loads and generators
Ygen_cut = 1./(1j*Xd1(iG_cut)); Bgen_cut = imag(Ygen_cut); % cuttable generator admittance
Yld_cut = (bus(iLD,4)-1j*bus(iLD,5))./(bus(iLD,8).^2);     % cuttable load admittance
Gld_cut = real(Yld_cut); Bld_cut = imag(Yld_cut);
Pgen_cut = gen(iG_cut, 2); 
Pld_cut = bus(iLD, 4);                                    % Pgen_cut & Pload_cut 
%% Initial state of generators
% No salient pole effect: Xq = Xd'
V_gen = bus(iG,8).*exp(1j* bus(iG,9));
% I(0) = I_x(0) + j I_y(0) = conj(S(0)) / conj(V(0));
% E_Q(0) = E_Qx(0) + j E_Qy(0) = V(0) + (Ra + j Xq)I(0).
% {x, y} => rotate (delta-pi/2) i.e.[sin,-cos; cos,sin] => {d, q}
% E'q = Vq(0) + Ra*Iq(0) + X'd*Id(0); Pm(0) = Pe(0) = P(0) + |I(0)|^2 *Ra.
I_gen = (gen(:,2) - 1j* gen(:,3)) ./ conj(V_gen);
E_Q = V_gen + 1j* Xd1 .* I_gen; 
delta = angle(E_Q); omega = ones(nG, 1);
V_dq = V_gen .*exp(1j* (pi/2-delta)); I_dq = I_gen .*exp(1j* (pi/2-delta));
E1_q = imag(V_dq) + Xd1.*real(I_dq);                       % real vector
Pm = gen(:,2);
% auxiliary variables
EXdTJ = E1_q./(Xd1.*TJ);
EXd = E1_q./Xd1;

end